#####################################################################
# Outcomes
#####################################################################

rm(list=ls())

sink("~/Dropbox/Geographic Natural Experiments/02_analyses/06_match_outcomes.txt")

library(Hmisc)
library(exactRankTests)
library(ggplot2)
library(reshape2)
library(xtable)
library(stargazer)
library(sensitivitymw)
library(sensitivitymv)
library("ggplot2")
library(gridExtra)

##############################################
# Read data and prepare basic output
##############################################

load("~/Dropbox/Geographic Natural Experiments/02_analyses/01_match_1000m.RData")
load("~/Dropbox/Geographic Natural Experiments/02_analyses/02_match_city_1000m.RData")
load("~/Dropbox/Geographic Natural Experiments/02_analyses/03_match_state_1000m.RData")
load("~/Dropbox/Geographic Natural Experiments/02_analyses/04_match_country_1000m.RData")
names(d)

##############################################
# Results perceptions
##############################################

# generate datasets for test
test_d_match = data.frame(d_match$turnout08[d_match$t_ind==1],d_match$turnout08[d_match$t_ind==0])
colnames(test_d_match) = c("treated","control")
test_d_match_city = data.frame(d_match_city$turnout08[d_match_city$t_ind==1],d_match_city$turnout08[d_match_city$t_ind==0])
colnames(test_d_match_city) = c("treated","control")
test_d_match_state = data.frame(d_match_state$turnout08[d_match_state$t_ind==1],d_match_state$turnout08[d_match_state$t_ind==0])
colnames(test_d_match_state) = c("treated","control")
test_d_match_country = data.frame(d_match_country$turnout08[d_match_country$t_ind==1],d_match_country$turnout08[d_match_country$t_ind==0])
colnames(test_d_match_country) = c("treated","control")

# point estimate 
pe_regular_gamma1 = as.numeric(sprintf("%.5f",senmwCI(test_d_match,gamma=1,method="t")$PointEstimate))[1]
pe_city_gamma1 = as.numeric(sprintf("%.5f",senmwCI(test_d_match_city,gamma=1,method="t")$PointEstimate))[1]
pe_state_gamma1 = as.numeric(sprintf("%.5f",senmwCI(test_d_match_state,gamma=1,method="t")$PointEstimate))[1]
pe_country_gamma1 = as.numeric(sprintf("%.5f",senmwCI(test_d_match_country,gamma=1,method="t")$PointEstimate))[1]
pe_regular_gamma1
pe_city_gamma1
pe_state_gamma1
pe_country_gamma1

# p-values gamma 1
pval_regular_gamma1 = as.numeric(sprintf("%.5f",senmw(test_d_match,gamma=1,method="t")$pval))
pval_city_gamma1 = as.numeric(sprintf("%.5f",senmw(test_d_match_city,gamma=1,method="t")$pval))
pval_state_gamma1 = as.numeric(sprintf("%.5f",senmw(test_d_match_state,gamma=1,method="t")$pval))
pval_country_gamma1 = as.numeric(sprintf("%.5f",senmw(test_d_match_country,gamma=1,method="t")$pval))

# p-values gamma 1.1
pval_regular_gamma1.1 = as.numeric(sprintf("%.5f",senmw(test_d_match,gamma=1.1,method="t")$pval))
pval_city_gamma1.1 = as.numeric(sprintf("%.5f",senmw(test_d_match_city,gamma=1.1,method="t")$pval))
pval_state_gamma1.1 = as.numeric(sprintf("%.5f",senmw(test_d_match_state,gamma=1.1,method="t")$pval))
pval_country_gamma1.1 = as.numeric(sprintf("%.5f",senmw(test_d_match_country,gamma=1.1,method="t")$pval))
pval_regular_gamma1.1
pval_city_gamma1.1
pval_state_gamma1.1
pval_country_gamma1.1

# p-values gamma 1.2
pval_regular_gamma1.2 = as.numeric(sprintf("%.5f",senmw(test_d_match,gamma=1.2,method="t")$pval))
pval_city_gamma1.2 = as.numeric(sprintf("%.5f",senmw(test_d_match_city,gamma=1.2,method="t")$pval))
pval_state_gamma1.2 = as.numeric(sprintf("%.5f",senmw(test_d_match_state,gamma=1.2,method="t")$pval))
pval_country_gamma1.2 = as.numeric(sprintf("%.5f",senmw(test_d_match_country,gamma=1.2,method="t")$pval))
pval_regular_gamma1.2
pval_city_gamma1.2
pval_state_gamma1.2
pval_country_gamma1.2

matching = c("regular matching","template:city","template:state","template:country")
pe = c(pe_regular_gamma1,pe_city_gamma1,pe_state_gamma1,pe_country_gamma1)
pval_gamma1 = c(pval_regular_gamma1,pval_city_gamma1,pval_state_gamma1,pval_country_gamma1)
pval_gamma1.1 = c(pval_regular_gamma1.1,pval_city_gamma1.1,pval_state_gamma1.1,pval_country_gamma1.1)
pval_gamma1.2 = c(pval_regular_gamma1.2,pval_city_gamma1.2,pval_state_gamma1.2,pval_country_gamma1.2)
table = data.frame(matching,pe,pval_gamma1,pval_gamma1.1,pval_gamma1.2)
table

stargazer(table, 
          type = "text", 
          colnames = FALSE, 
          summary = FALSE, 
          title="Permutational t-test and sensitivity analysis",
          digits=2, 
          rownames=FALSE, 
          float = TRUE, 
          float.env = "table", 
          table.placement = "H", 
          #column.sep.width = "-30pt",
          out="~/Dropbox/Geographic Natural Experiments/03_manuscript/tables/t4_results.tex")

